arXiv:1506.02856vl [cond-mat.mes-hall] 9Jun2015 
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We theoretically examine the role of Kerr nonlinearities for graphene plasmonics in nanostructures, specifically 
in nanoribbons. The nonlinear Kerr interaction is included semiclassically in the intraband approximation. The 
resulting electromagnetic problem is solved numerically by self-consistent iteration with linear steps using a real- 
space discretization. We derive a simple approximation for the resonance shifts in general graphene nanostructures, 
and obtain excellent agreement with numerics for moderately high field strengths. Near plasmonic resonances 
the nonlinearities are strongly enhanced due to field enhancement, and the total nonlinearity is significantly 
affected by the field inhomogeneity of the plasmonic excitation. Finally, we discuss the emergence of a plasmonic 
bistability which exists for frequencies redshifted relative to the linear resonance. Our results offer new insights 
into the role of nonlinear interaction in nanostructured graphene and paves the way for experimental investigation. 


Nonlinear optical effects, 1 - 2 facilitated by strong light-matter 
interaction, are indispensable in modern photonics. Indeed, a 
host of phenomena and applications arise at sufficiently high 
field-strengths, owing to superlinear photon-photon response 
mediated by strong light-matter interaction, ranging from fre¬ 
quency conversion, through all-optical phase-modulation, to 
ultra-fast switching, and is pursued in a broad range of plat¬ 
forms. 3-5 

A perennial challenge in the discipline is to achieve signif¬ 
icant nonlinear interaction at ever smaller excitation powers 
and interaction volumes, whilst maintaining in-situ tunability 
and control. In achieving this goal, the field of plasmonics, 
describing the strong hybridization of the free electromagnetic 
field with collective oscillations of conduction electrons, sug¬ 
gests several promising avenues. 6 In particular, the extreme 
local field enhancements inherent to plasmonic excitations 
amplify intrinsic nonlinearities considerably, allowing large 
effective nonlinearities. 

Nevertheless, plasmonic field-enhancement is fundamen¬ 
tally limited by intrinsic Ohmic losses even in noble metals. 
The advent of the two-dimensional material graphene has gar¬ 
nered significant interest in the plasmonic community, 7 - 3 in 
part due to extremely large electron mobilities 9-11 and con¬ 
comitant extraordinary plasmonic field-enhancements, 12 ex¬ 
ceeding even the very large enhancements known from metal- 
plasmonics. Furthermore, graphene has attracted much inter¬ 
est also for its exceptional intrinsic nonlinear properties both 
theoretically 1 ’ -l6 and experimentally. 17-19 Building on this 
compound-fortuity, a body of research is rapidly emerging at 
the crossroad of nonlinear plasmonics and graphene. 20-26 

Very recently, the role of Kerr nonlinearities in infinitely 
extended graphene has been studied, 24 notably establishing 
the existence of bistable solutions. In this paper we study theo¬ 
retically an analogous Kerr nonlinearity but in nanostructured 
graphene, specifically in nanoribbons - wherein plasmons, 
unlike in the extended counterpart, are readily excited with¬ 
out momentum-matching concerns, e.g. by normally incident 
plane waves. We report an induced nonlinearity which is 
strongly affected by the degree of inhomogeneity of the elec¬ 
tric fields of the plasmon - a feature which is absent in the 
corresponding extended system. Furthermore, we derive a 
simple perturbative expression for the nonlinear resonance 
shifts in general graphene nanostmctures, and show that it is in 


excellent agreement with full self-consistent calculations for 
moderately high field strengths. Finally, we discuss the emer¬ 
gence of plasmonic bistability in nanoribbons under normally 
incident plane-wave excitation. First, however, we introduce 
the two basic components needed for a nonlinear treatment of 
graphene nanostructures, namely a material response model 
and an exposition of the resulting electromagnetic problem. 

Material response. For photon energies hoj low compa¬ 
rable with the Fermi energy e F , the response of graphene is 
reasonably approximated by neglecting interband transitions. 
In this case, the intraband response can be derived from the 
Boltzmann equation. To third order in the perturbing field the 
Kerr-corrected conductivity, i.e. the response oscillating at the 
perturbing frequency, is 24 


c(r) = cr (1) [l - |E(r)| 2 /£ ( 2 3i ], (1) 

expressed in terms of the linear intraband conductivity cr (1) = 
ie 2 e F /7 Th 2 (aj + iy) with loss-rate y, and a third-order charac¬ 
teristic field £ 2 3) = (8nr 2 ) )/(9« 2 )£’ 2 at linearly related to the 
saturation field E sM = e, oj/ev, through a loss-modified fre¬ 
quency ttr 2 3| = (o) + -fiyiioj - iy). Since the Kerr correction 
is of the self-focusing type, 2 - its usage in finite structures 
with inhomogeneous fields must be augmented to include a 
saturating mechanism, or else suffer nonphysical run-away 
self-focusing. 2 Here we adopt the well-known two-level satu¬ 
ration model, or, in other words, the [0/2] Pade approximant 
of <r(r) consistent with Eq. (1) 


cr(r) 


PinCr) 

1 + |E(r)| 2 /E 2 


+ cr ( 3 J2r (r). 


( 2 ) 


This model reproduces the third-order result of Eq. (1) in the 
|E(r)|/£ sat <s 1 limit, while crucially exhibiting a sensible 
behavior beyond this limit as well. 3 Lastly, we include in 
Eq. (2) a term <x< 3)2 ,,(r) to account for a high-field loss mecha¬ 
nism through two-photon absorption via the phenomenological 
prescription suggested by Gorbach, 23 via the dissipative correc¬ 
tion cr (3)2y (r) = -io' 27 cr (ll |E(r)| 2 /£' 2 at with a 2y * 0.1 estimated 
from measurements. 19 

Before proceeding we briefly discuss the limitations of the 
material response assumed in Eq. (2). Firstly, the disregard 
of interband effects limits our consideration to energies suf¬ 
ficiently below ~2e, . Secondly, nonlocality, 29 edge-states, 30 
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and more generally atomistic features 25,26,31,32 are excluded, 
though they are important at small feature sizes. Consequently, 
we restrict our considerations to nanostructures of character¬ 
istic dimensions > 25 nm where these effects only weakly 
perturb the intraband approximation. 

Interacting response. In the quasistatic limit, the interacting 
response of graphene can be deduced from three elements; 
the Coulomb law, the continuity equation, and the current- 
field relationship as specified by a conductivity-model. For a 
nanostructure defined by a two-dimensional domain Q (e.g. at 
Z = 0), these elements combine to form an integro-differential 
equation for either the induced density or the total potential 
<p( r). Here we choose the latter 8 

0(r)=-—d 2 r' V(r,r')V'- [<x(r')V'</>(r')], (3) 

AnsaUiW J n 

expressed in dimensionless coordinates r (,) = [x ( ' , ,v ( ' 1 , z] 1 
normalized by a characteristic length W, with the Coulomb 
interaction V(r, r') = |r - r'| -1 , and with differential operators 
V' = [d X ',dy'Y. The conductivity cr(r) implicitly depends 
on frequency - and in a nonlinear treatment also on the total 
field E(r). The spatial dependence of the conductivity can 
be conveniently expressed via a dimensionless occupation 
function /( r) = cr(r)/(cr (1) ) with (cr (1) ) denoting the average 
linear conductivity across D. Introducing operators Vg(r) = 
f dr' V(r, r')g(r') and Dg(r') = V'- [/(r')V'g(r')] casts Eq. (3) 
as an eigenvalue problem for the composite operator VD 

A0(r) = VDcfi(r), (4) 

with eigenvalues A = 47r£o<uW/i(cr (1) ), dictating the permitted 
eigenfrequencies a>. Operators V and D find simple matrix- 
forms in a discretized real-space basis in both the general 2D 
case as well as in the ID ribbon case, see Supplemental Ma¬ 
terial (SM). The operator D is constructed so as to account 
explicitly for a boundary condition of vanishing normal cur¬ 
rent (or, equivalently, for the conductivity-discontinuity) at 
the boundary <90. In the presence of an external potential 
0 ext , the eigenvalue problem in Eq. (4) becomes an inhomoge¬ 
neous equation through the addition to the right-hand-side of 
a source-term A<f> ext ( r). To solve the nonlinear problem, with 
<x(r), and hence /( r) and D, depending on the total electric 
field locally, we solve the nonlinear system iteratively until 
self-consistency is reached, exploiting at each iteration-step 
the computational efficiency associated with linear systems, 27 
see SM. 

With the formal premise established, we next specialize to 
the case of nanoribbons, translational ly invariant along y and 
of finite extent W along x; a system which has already attracted 
much attention in the linear case. 29,32-34 As a consequence of 
translational symmetry, eigensolutions can be expanded in a 
momentum basis according to 0(r) = </>(x, z ) exp(i^ny). Of key 
interest is the evolution of the eigenenergies with momentum 
k t (here dimensionless; conventional units via kJW), i.e. the 
dispersion relation Ha> n (k tl ) - and subsequently the response of 
the system to external fields. 

Eigenmod.es and nonlinear dispersion. For low field 
strengths, i.e. in the linear regime with /(r) independent of 
E(r), the eigenmodes A„(k 3 ) of Eq. (4) are solely geometry 



FIG. 1 (color online), (a) Dispersion relation of a single nanoribbon. 
Ribbon-averaged field strength (|E(r)|) ranges from negligible (black), 
i.e. linear, through lx 10 5 V/cm to 4x 10 5 V/cm (lightest blue) in steps 
of 0.5 x 10 5 V/cm (increasing along arrow). For the first hve (|E(r)|), 
we indicate in dashed red the corresponding analytical estimate, see 
Eq. (5). For the monopole, only the linear calculation is shown. The 
region of significant interband-modification is illustrated in shaded 
gray. Inset schematically depicts a single graphene nanoribbon, (b) 
Field intensity, |E(r)|, contour maps for the case (|E(r)|) = 4 x 
10 5 V/cm and k t = 0. Colormap ranges from maximal (dark) to 
minimal (light) logarithmically, with contours separated by factors of 
1.5, 1.75, 2, and 2.25 for dipole, tripole, quadrupole, and pentapole 
cases, respectively. Sparklines below maps depict the variation of 
|/(r)| along the ribbon, with maximal and minimal values indicated. 


dependent - but scale invariant - with associated eigenfrequen¬ 
cies «„(£„) dictated by A n (k 3 ) = 47reo<u„(^!i)W/i(cr a) ), allowing 
in the linear intraband approximation the simple scaling re¬ 
lation «„(£„) ^ (2^) _1 sj-A„{k^e 2 e ? ls^W . 29,33 Under signifi¬ 
cant nonlinear interaction, however, the eigenvalues A n (k 3 ) are 
field-dependent, and, by extension, scale-dependent due to the 
self-consistent nature of the problem. In Fig. 1(a) we investi¬ 
gate the dispersion relation of the first few eigenmodes of a 
single W = 50 nm nanoribbon for different ribbon-averaged 
field strengths (|E(r)|) = f n dx |E(x)|. The most apparent 
impact of nonlinearity is a redshift of all resonances. This 
is readily appreciated from the negativity of the Kerr correc¬ 
tion. Indeed, the shift can be well-approximated by perturba¬ 
tion theory for any general structure: denoting by luo',,' and 
Ef the linear response eigenenergies and eigenfields [with 
(|E™ ) (r)l) = (|E(r)|)] the nonlinear eigenenergies are, to lowest 
order, approximately (see SM) 


(5) 


with the averages taken over r e O. The approximation is 
excellent for moderately high fields, see dashed red lines of 
Fig. 1(a), though, naturally, inaccurate for the largest con- 
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sidered fields due to the disregard of the self-consistent as¬ 
pects of the nonlinear perturbation. The approximation also 
reveals the important role played by the inhomogeneity of 
E(r), or equivalently /( r), for the nonlinear strength since 
<|E <0) (r)[ 4 ) + (|E IO) (r)| 2 ) 2 for inhomogeneous fields. 

In Fig. 1(b) we explore this point further, by depicting 
the modal character and inhomogeneous nature of the plas- 
monic modes. The modal labels are chosen from the per¬ 
spective of the induced charge density, p(x), of the mh mode, 
with the monopole, dipole, tripole, quadrupole, and pentapole 
(n = 0,1,2,3, and 4, respectively) exhibiting n nodes of p(x). 
Modes of even n are optically dark, owing to a vanishing 
dipole moment, and remain optically dark also under nonlin¬ 
ear perturbations (which preserves the system symmetry). The 
monopole violates charge conservation along x [but not along 
(x, y) for k, + 0], is optically dark, and consistently does not 
converge at higher fields; as a consequence, we depict only 
its linear dispersion. The variation of the occupation function 
/(r) under large fields is highlighted in the insets of Fig. 1(b). 
The strong spatial variation of the conductive profile, up to 
50% for the considered (|E(r)|), is a direct consequence of the 
strongly inhomogeneous nature of plasmons. Despite these 
significant spatial variations of /(r), the corresponding modal 
profiles away from the ribbon are highly similar under linear 
and nonlinear circumstances, since their character is dictated 
chiefly by the nodal character of p(x). 


I/V = 25 nm 


tV = 50 nm 



FIG. 2 (color online). Field enhancement (|E(r)|)/£o as a function 
of energy Tito, for varying incident field strengths E 0 (as indicated 
above each spectrum). Each spectrum is offset vertically by 5 units. 
Two ribbon widths W = 25 nm and 50 nm are examined. Regions of 
bistability are delimited by dashed arrows which indicate the ramping 
direction. Material parameters are as in Fig. 1. 


Plane-wave excitation and bistability. Having considered 


the dispersion of eigenmodes, we next turn our attention 
to the response of the system due to a normally incident 
plane wave, polarized along x, i.e. E ex t(z = 0) = £yx and 
<pext(z = 0) = -EqxW, corresponding to vanishing k ]v In addi¬ 
tion to the power absorbed from the incident wave, the induced 
and total electric fields are of primary interest - here we focus 
on the latter. For reasons of numerical efficiency, and as we 
shall see, physical necessity, we compute for each separate en¬ 
ergy the response by an initial linear calculation, followed by 
a ramping of the incident field strength. Specifically, for fixed 
tiu), we consider a ramp-array {Eo t „}^ =l with Eo.h+i > £o,« and 
with £0 [ sufficiently small to be considered a linear perturba¬ 
tion. Starting from £ 0,1 we compute associated solutions and 
proceed, generally, to field strength Eo.n+i with initial guesses 
on / and tf> obtained from the nth solution. This defines the 
upward ramp, corresponding to slowly turning the incident 
intensity up. Upon reaching n - A' we invert the procedure 
and follow a downward ramp, in the pattern £o,„ —» Eo.n-i, 
corresponding to slowly turning the intensity down. 

In Fig. 2 we examine the spectral response of ribbons of 
widths W — 25 nm and 50 nm under different excitation 
strengths, i.e. under varying £ 0 . For moderately high Eq 
the linear Lorentzian resonance is asymmetrically perturbed, 
slightly broadened, and redshifted. Furthermore, the upward 
and downward ramps to Eq give identical spectra. As Eq is 
increased further, these perturbations intensify. However, in 
certain frequency ranges the response to upward and down¬ 
ward ramps toward £0 differ (regions delimited by dashed 
arrows); a trademark of bistability. Similar features were dis¬ 
cussed for extended graphene in Ref. 24 using the Kerr model 
of Eq. (1) and in Ref. 25 for finite systems using a phenomeno¬ 
logical anharmonic model. The primary extension here is the 
full self-consistent accounting of the inhomogeneous nonlinear 
conductive profile arising in nanostructured systems. For com¬ 
parison, we note that the bistability studied here is achievable 
at much larger energies than in the extended system, where it 
is restricted to Tito < V4/3ct' fs e F under normal incidence (with 
a fs = e 2 /AnEohc). 1 ' Here, bistability is evident in the dipole 
mode for both W — 25 nm and 50 nm, but also visible for 
the quadrupole mode for W = 50 nm. In both cases, the area 
traced by the bistable region initially increases with Eq and 
then decreases due to saturation and increased absorption. 

The history dependence of the response is further examined 
in Fig. 3(b), depicting hysteresis curves of Eq vs. (|E(r)[) 
(normalized to the frequency-dependent saturation field) at a 
selection of fixed frequencies as indicated in the linear spec¬ 
trum of Fig. 3(a). At energies far from the linear resonance 
at ho/"’ the response (|E(r)|) relates linearly with £ 0 . As the 
energy is increased towards ho/", a nonlinear discrepancy de¬ 
velops with increasing Eq which eventually gives way to a 
discontinuous jump at a critical field strength £j, indicated 
for a selected energy in Fig. 3(b). As Eq is reduced on the 
downward ramp, its response initially traces out that of the 
upward ramp, but departs from its upward correspondent after 
£ () and eventually undergoes a discontinuous jump at £ () after 
which the initial path is retraced. The hysteresis area, indi¬ 
cated by shaded areas, increases with positive to m - to (though 
El similarly increases, delaying the onset of hysteresis), but 
vanishes for to > to <0> due to the redshifting of the resonance 
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FIG. 3 (color online). Hysteresis arising from bistable behavior in 
a W = 25 nrn nanoribbon excited by a plane wave E 0 x (material 
parameters are as in Fig. 1). (a) Linear response field-enhancement 
spectrum versus energy. Selected energies are highlighted by colored 
markers, and the linear resonance energy Tia> m is labeled explic¬ 
itly. (b) Hysteresis curves at fixed energy [corresponding colorwise 
to those highlighted in (a)] for total field (|E(r)|) versus incident 
field Eq. Bistable regions are indicated by shading and delimited by 
energy-dependent low- and high-point field strengths E^. (c) Inten¬ 
sity maps of the induced electric field Re[ E'" d (x, z) |. Colorscale is 
identical across the four maps, ranging from positive (red), through 
zero (white), to negative (blue) in a symmetric range. Absolute mag¬ 
nitudes are scaled logarithmically for intelligibility. Frame color 
indicates association with energies in (a). Field strengths in the high- 
field maps are specified by corresponding triangles in (b). Sparklines, 
defined as in Fig. 1(b), indicate the range and variation of |/(r)|. 


with Eq. The onset of bistability is reached for incident field 
strengths considerably below E sat ; this fortuity is of course a 
result of plasmonic field-enhancement of the total field. 

Lastly, we comment on the field profiles of the excitations. 
First we highlight the linear response at energies just below 
and above ha> t0) , indicated by red and green markers in Fig. 3(a) 
and 3(c). The field profile exhibits a well-known n phase-shift 
between the two energies, a result which can be appreciated 
e.g. by inspection of the linear harmonic-oscillator polarizabil¬ 
ity a(u>) oc [(o/ 0) ) 2 - cj(cl) + i~y)] “ 1 which exhibits a sign-change 
of its real part as a> traverses a> m : as a result, the induced dipole 
p(o )) = a{u))EQ changes sign for u> S a> <0) , and correspondingly 
so for the induced fields. A similar sign change is observed in 
the bistable comparison, see black-framed modes in Fig. 3(c). 
Again, the origin of the sign change can be appreciated from a 
polarizability consideration by including a third-order anhar- 
monic term to the harmonic oscillator model, 2 '’ see SM. 


Summary and discussion. In this paper we have analyzed 
the impact of Kerr nonlinearity on the plasmonic response of 
graphene nanostructures, specifically for nanoribbons. The 
key distinction of nanostructures compared to the correspond¬ 
ing extended system arises from the strongly inhomogeneous 
fields of localized plasmonic excitations, which in turn incur 
an inhomogeneous conductive profile. We have derived a 
simple analytic expression, Eq. (5), which approximates the 
nonlinear resonance shifts, while accounting for both inhomo¬ 
geneity and overall amplitude of the nonlinear perturbation. 
The characteristic field of the Kerr nonlinearity in graphene is 
the saturation field E sat . However, significant nonlinear interac¬ 
tion can be achieved near plasmonic resonances even for much 
weaker incident fields owing to plasmonic field enhancement. 
Finally, we discussed the existence of a plasmonic bistability 
in nanoribbons under normal incidence. 

The applications of optical bistabilities are well-known 
and long-pursued, 1,2 with implications particularly in optical 
switching. Indeed, a range of platforms have been scrutinized 
for this purpose, in recent years e.g. in photonic crystal cavities 
(PCC) where nonlinearities are enhanced by large 2-factors 
and light-slowdown. 3 Whether graphene can further the state- 
of-the-art in this mature field remains to be seen. 35 We expect, 
however, that a very profitable avenue for progress exists in 
hybrid approaches, utilizing e.g. PCC and graphene in unison 
- as has in fact been explored experimentally, 19 albeit with¬ 
out taking advantage of the resonant plasmonic nonlinearity 
described herein. A simultaneous tuning of both cavity and 
plasmonic resonance should allow for maximal nonlinear man¬ 
ifestation in such systems. Advances in this direction requires 
improved understanding of nonlinearities in nanostructures; 
the present work constitutes one such effort. Several features, 
however, remain unexplored, underscoring the fertility and 
richness of the field. For example, from a semiclassical per¬ 
spective, barring atomistic approaches, 23,211 questions remain 
relating to the role of interband nonlinearities, 16 nonlocality, 
and the effective role of edge states. 

In closing, we mention a final question of singular practical 
relevance, namely damage thresholds. So far, to the best of 
our knowledge, measurements do not exist in the infrared, but 
in the optical domain 36 38 the reported thresholds fall in the 
rather broad range from ~ 10 6 V/cm in fs-pulsed operation 
to just ~ 10 4 V/cm for hour-long continuous wave operation. 16 
For comparison, the saturation field at tia> — e T — 0.2 eV is 
£ s at ~ 6.7 x 10 5 V/cm. Though direct comparison is impos¬ 
sible, in part due to frequency range, pulse conditions, and 
the uncertain impact of field enhancement, this highlights that 
even resonantly enhanced nonlinearities in graphene walk a 
narrow road - not unlike previous contenders for large nonlin¬ 
earities. Given the promising results presented herein, how¬ 
ever, we believe the journey will be worth the effort. 
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SUPPLEMENTAL MATERIAL 


I. ITERATIVE PROCEDURE FOR NONLINEAR PROBLEM 

We here discuss an iterative approach to solving the nonlinear equation 

Hi r) = #ext(r) + VD[/[0]]0(r), (SI) 

which is essentially just the driven correspondent of Eq. (4), and where we have emphasized the dependence 
of D on /(r) through /(r). The problem is evidently nonlinear, but can be solved efficiently by iteration 
with only linear algebra at each step. We follow the usual iteration scheme, as e.g. also used previously in 
the studies of bistability in dielectric waveguides: si 

1. Compute a linear solution based on an initial guess of / = fi m , i.e. solve Eq. (SI) with D [/[/]] —> 
D[/ = /in;]. Denote the obtained solution as (jp\ Set the iteration step m = 0. 

2. Calculate the mth guess at the occupation function / [ml from the potential 

3. Compute the (m + l)th iteration by solving the linear system = A(p ext + VD[/ [ml ]/ [m+1] . 

4. Iterate steps 2 and 3 until convergence, otherwise update iteration step m —> m + 1. 

We impose convergence criteria corresponding to the simultaneous fulfillment of (with tol = ICC 5 ) 

max U [m+1] (r) - / [ " ,] (r)|/max U [ffl] (r)| < tol, (S2a) 

rsn rEfi 

max |/ [m+1] (r) - / [m] (r)|/max |/ [ffl] (r)| < tol, (S2b) 

being of standard type for iterative approaches to nonlinearity. 1 In all considered cases the iterative 
procedure converged after at most several hundred iterations. One exception should be mentioned however; 
the dipolar eigenmodes at field strengths 3 x 10 5 V/cm and 3.5 x 10 5 V/cm failed to converge after 1250 
iterations for k ]t > 5 and are consequently absent in Fig. 1 for these momenta. This could likely be 
remedied by a more elaborate stepping procedure, though such investigations have not been pursued 
further in this work. 

Two additional extensions of the simple iterative scheme described above are employed. Firstly, for 
numerical stability we apply a linear mixing scheme for updating guesses on /, specifically we use 
D[/«] with /[” ] = (1 - x)/ [m ~ 1] + £nix/ [m] in Step 2 (mixing parameter ^ = 0.275) rather than 

the unmixed D[/ [m1 ]. Secondly, the initial guess /,„ is always taken from the previous field strength in 
ramping scenarios. This provides a significant numerical speed-up and, crucially, allows us to investigate 
hysteresis and bistability. The initial guess at the first field strength is naturally / n ; = 1. 

For eigenmodal calculations where </> ext = 0, we normalize (/>„ at each iteration to impose the desired 
ribbon-averaged field strength (|E(r)|), and in addition determine u>„ from A„(a>„) by numerically solving 
the equation in the complex frequency-plane. 


II. MATRIX REPRESENTATION OF V AND DINA DISCRETIZED BASIS 

We here elaborate the reduction of the differential and integral operators D and V to matrix representations 
D and V using an equidistant discrete basis. Specifically, we discuss the ID ribbon case, although the 
generalization to general 2D restrictions is straightforward. Specifically, we imagine a system in the 
xy-plane, translationally invariant along y and with finite extent along x. For simplicity, we assume just 
a single ribbon, such that x is limited to the simple domain a: e [0,1]. Furthermore, as the operators 
necessarily act on a potential /(r), we impose translational invariance along y by the decomposition 
/(r) = <p(x )e tfc|,y . 

Starting with the differential operator D, we consider its operation onto <p(r), which takes the form 
D/(r) = d x [f(x)d x tp(x)]e lhiy - f(x)(p(x)e' k,y . By extension, we define the operation of D onto the 
single-variable function 4>(x) through D/(x) = d x [f(x)d x (t>(x)] - k 2 f(x)(b(x). To proceed, we introduce 
a discretization of the x-coordinates as with associated values <pj = (Mx/) and / = f(xj) (we take 
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FIG. SI (color online). Sketch of the discretization approach applied to a single ribbon. 


N - 150, being well-converged in all considered cases). Though not strictly necessary, we assume 
equidistant Xj with constant spacing Xj+\ - Xj = a, see Fig. SI. 

The matrix elements Dj\ of the finite-element representation of D is then defined by D<pj = EzTJ/z^z- The 
elements can be deduced using finite differences at the midpoints. Specifically, using central differences 
dJifMj] ~ a ~ 1 ( m / _ m j- 1) where mj defines midpoint-values of the function mix) = f{x)d x <p{x) such 
that mj — {2a)~ l {fj + \ + fj\<f>j+ i - </>;), see Fig. SI. For all interior points, j e [2,N - 1], this then allows a 
decomposition of Dji as the tridiagonal matrix 

Dji = 2?[^7-u(^-i + fj) - tijiifj -1 + 2 fj + fj+ 1) + dj+iAfj + fj+i)] - fyrffj- (S3a) 

At the end-points j = 1 and j - N we explicitly account for boundary conditions. Specifically, we ensure 
a vanishing of normal current, equivalent to the condition d x <p(x) = 0 for x = 0 and x = 1. In turn, this 
forces mo = m N = 0, allowing 


Du = ^(/i + m-S u + 8 2 ,i) - 8 u tff u (S3b) 

Dm = 2^2 (/az-i + f N )(8 N -v - 8 n,i ) _ 8 N jk~f N . (S3c) 

As an alternative to taking explicit account of the boundary condition, one can allow a slightly larger 
x-range, and explicitly include points with /(r) = 0 outside r e Q - the step in /(r) at r e dCl then mimics 
an edge charge and accounts numerically for the boundary condition; such a procedure may be preferable 
in finite structures without any geometric symmetries compatible with a square grid. 

The integral operator V is similarly amenable to explicit expression on the equidistant grid. Specifically, 
letting V operate on a function g(r) = g(x)e l, |V one finds S2,S3 

Vg(r) = ■ dx' 2 K 0 (k„\x - x'\)g(x'), (S4) 

where k tl > 0 is assumed and with Kq denoting the zeroth order modified Bessel function of the second 
kind. Assuming a slowly varying g(x) and an equidistant {xj) then allows a matrix decomposition of V via 
Vg, = Ez Vjigj where S4 

r»Xi+a/2 

Vji = 2 I dx'K 0 (kJx j -x'\) = nY J x{K 0 (k ]] \x\)[L 1 (k l jx\)+l] + K 1 (k l \x\)L 0 (k l \x\)\, (S5) 

Jx,-a/2 x=Xj!±a/2 

with xji = xj - x; and Lq i denoting modihed Struve functions of zeroth and first order. 

A final detail which should be discussed is the special case k n = 0, where the kernel K f )(k [x - x'|) 
in Eq. (S5) diverges. Despite this divergence, finite and meaningful matrix elements can be retrieved 
by invoking charge conservation. Specifically, we note the small argument expansion A'q^Jx - x'[) ~ 
- ln(|x-x , |)-ln(£||)+a' where a = ln(2)-y EM (y EM is the Euler-Mascheroni constant). S4 The x'-independent 
term - lnf/r,) + a gives a contribution [- ln(^,|) + a] f dx' g(x') to Eq. (S4) and appears divergent as k tl —> 0. 
Nevertheless, this contribution vanishes for the functions g(r') of relevance since they always represent 
induced charges [as evident from Eq. (3)] and obey charge conservation f dx' g(x') - 0. As such, the 
k ]t = 0 case can be calculated by simply letting Ko(k t \x - x'|) —> -ln(|x - x'|) in Eq. (S5), s 1 yielding 
Vji = -2 Yj S =± s(xji + sf) Infix;/ + sf |) for k tl = 0. 

This concludes the real-space discretization approach for reduction of the abstract operator equation of 
Eq. (4) into a matrix equation A<p - \ZD</> with <f> denoting the vector form of cf> r 
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III. PERTURBATION ESTIMATE OF THE NONLINEAR SHIFT OF EIGENFREQUENCIES 


We here provide the derivations that allow the approximate result of Eq. (5). As we explain below, the 
approach relies on the formulation of a Hermitian eigenproblem followed by application of standard 
perturbation theory to a spatially inhomogeneous problem. 

The compound operator VD defined in Eq. (4) is - though numerically practical - inconvenient for 
analytical considerations, because it is not symmetric. However, the problem can (of course) be cast 
as a Hermitian eigenproblem with eigenvalues A„ [though, strictly speaking, only for real, positive 
occupation functions /(r), which we restrict our analysis to here], as also noted recently in Refs. S5 
and S6. Specifically, consider the application of the scaled gradient operation - y/( r)V onto Eq. (3): 

- A V7(r)V0(r) = V7WV £d 2 r' V(r, r')V' ■ { y/f(?) [ - VflFjV'tfr')]}. (S6) 


Defining the scaled in-plane field £(r) = - \jf(r)V6(r ) and manipulating further allows 

Ag(r)= V/WV f d 2 r'V(r,r')V'-[Jf(Pjf(r')] 

Jn 

= VTOOvj £d 2 r' V' • [V(r, r') V/TOCO] - £d 2 r' [V'V(r, r')] • [ r')] 

= - V7^)V £d 2 r' V/(P)[V'y(r, r')] • £(r') 

= - f dV V/(r)/(r')[V®V'y(r,r')]^(r') 

J n 

with associated steps a - c explicated below for convenience: 


(S7) 


a. Application of chain rule to expand integrand. 

b. The first integral term in step a vanishes, as can be deduced by application of the divergence theorem 
which transforms the term to <fc n V(r, r') \jf(v')\^(r') ■ n']. The integrand vanishes for all r' e <9Q 
due to the no-spill boundary condition on the induced current which forces £(r') • n' = 0 on r' e 8L 2. 

c. The term Jf( r)V is taken under the integral sign. V operates on r and hence only on V(r, r')- The 
operation V{[V'V(r, r')] • v(r')} is rewritten in the equivalent outer-product form [V® V'V(r, r')]v(r') 
with elements [V ® V'],-,- = d, d/. 

J i j 

We then define the operator M by its action on a field-ket |f) [where, as usual, (r|£) = £(r)] 

<r|M|f)= f d 2 r' V/(r)/(r')[V ® V'V(r, r')^(r'), (S8) 

Jn 

with associated eigenspectrum (-A n ) and \^ n ): 


(-A n M n ) - M|£,>. (S9) 

The operator M is evidently symmetric, positive semi-definite, and thus Hermitian. Aaccordingly, the 
eigenspectrum {-/!„) is non-negative and real; and the eigenkets \g„) are orthogonal = 6 n ri{£n\£n) 

and span the solution space for r e £1 

With these facts established, we can now discuss a perturbation treatment. Specifically, we consider the 
simple case where /( r) = f"" + ci/ ll, (r) for re!) with “groundstate” f° = 1 and perturbation / ' with 
strength 6. The corresponding expansion of M = M (0) + c)M !l> + 0(6 2 ) is found by expansion of Eq. (S8), 
yielding 


(r|M IO) |£> 



1 V'V(r, r')]^(r'), 

(1) (r) + / (1, (r')][V®V'y(r,r')]^(r'). 


(SlOa) 


(SlOb) 
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Since M is a Hermitian operator usual perturbation theory applies. S/ Specifically, for a “groundstate” 
eigenspectrum {-/C, the leading-order correction to the perturbed eigenvalue A„ = 4™ + 6A^' + 0(5 2 ) 
is derivable by application of Eqs. (S10) [by using the (r, r')-symmetry of the resulting equation] 


(ffiMiim _ m (C\no 

(C\C) " (C\C) 


For nonlinear purposes, we unfortunately do not know the exact perturbation / (1) as it should be determined 
self-consistently with the total field |£„). However, for low field-strengths this self-consistency can be 
neglected and we can approximate /[|£„)] — f\\£n)\ with |$° J ) referring to the electric field predicted by a 
linear calculation (at the desired field strength). For the Kerr-type nonlinearity of Eq. (1) the resulting 
correction is therefore [assuming vanishingly small loss and noting £"‘(r) = E <0) (r) for f <0> = 1] 


,o,9 /nd 2r IE<°>(r)| 4 

" ~~ " 8 £ 2 at J nd2r | E( „ )(r) | 2 


r> 9 (|E (0) (r)| 4 ) 


(Sll) 


with E sat similarly evaluated at the linear resonance frequency oj'n associated with d™. Finally, the result 
of the main text, Eq. (5), is obtained by invoking the relation between eigenvalues d„ and eigenfrequencies 
u> n together with the lossless intraband conductivity <x (1) (a;) ^ ie 2 e F /nti 2 u). 


IV. QUALITATIVE ANHARMONIC OSCILLATOR MODEL 

We review the basics of the simple anharmonic oscillator model, S8,S9 and discuss how it - in connection 
with a polarizability consideration - explains the n phase-shift observed for the bistable solutions in 
Fig. 3(c). 

In this qualitative model, we represent the induced dipole by a single (time-dependent) coordinate x, which 
obeys the simple equation of motion 

rax + myk = -efE^(t) - d x U(x), (S12) 

with an effective anharmonic restoring potential U(x) = \m(oJ 0) ) 2 x 2 - | max 4 , effective oscillator mass ra, 
linear resonance tu (0) , anharmonic parameter a (note that a > 0 in our case cf. sign of Kerr conductivity), 
and coupling factor /. We seek the solution that oscillates at e -1 "' in response to a perturbation EoCO = 
Eo(a>)e ~'‘ Jt , i.e. the Kerr response; we denote this term by x tl “ ,) (cj)e^ 10> '. Working with Eq. (S12) one finds 
(omitting declaration of ru-dependence) 

ra[(a/ 0 ’) 2 - «(w + iy) - 3a|x ,M | 2 |x' M = -efE 0 . (S13) 

The polarizability a m is linked to x (1 “ ) via the induced dipole = -ex <M = a^'Eo, allowing (ignoring 
loss, being nonessential for the present considerations) 

[(o/ 0) ) 2 - or - 3fle^ 2 [a ll “ , | 2 £^]a; (1 “ ) = e 2 f/m. (S14) 

For the bistable scenarios the term (w“ l! ) 2 - or is always positive, see e.g. Figs. 2 and 3. Depending 
on the magnitude of 3 ae^ 2 a 0cu) E^ relative to (of") 2 - or it is then clear that polarizability-solutions of 
opposing sign can arise, depending on the sign of the terms bracketed on the left-hand side of Eq. (S14). 
Furthermore, if we denote the positive and negative solutions a"l" 1 ' 1 and respectively, it can then be 
deduced by direct inspection of Eq. (S14) that \a ( + J> \ < |o£" , |. In other words, the induced dipole - and 
hence the induced fields - of the positive solution should be lower than its negative counterpart; upon 
identifying the lower branches of Fig. 3(b) with or ( + and vice versa for the upper branch, we see that this is 
exactly the case. As such, the anharmonic model describes not only the phase-shift, but also the magnitude 
interrelationship. Lastly, we mention for completeness that the anharmonic model describes also a third 
solution, which, however, is physically irrelevant as it is unstable (and correspondingly is not found in the 
iterative procedure employed in this study, nor in experimental investigation). 


* asger@mailaps.org 







10 


[51] X.-H. Wang, Finite Element Methods for Nonlinear Optical Waveguides , Advances in Nonlinear Optics, Vol. 2 
(Gordan and Breach Publishers, 1995). 

[52] I. Gradshteyn and I. Ryzhik, Table of Integrals, Series, and Products, 7th ed. (Elsevier Academic Press, 2007). 

[53] W. Wang and J.M. Kinaret, Phys. Rev. B 87, 195424 (2013). 

[54] National Institute of Standards and Technology, “Digital library of mathematical functions;' Release 1.0.9 of 
2014-08-29. 

[55] K.A. Velizhanin, Phys. Rev. B 91, 125429 (2015). 

[56] F.J. Garcia de Abajo and A. Manjavacas, Faraday Discuss. 178, 87 (2015). 

[57] S. Gasiorowicz, Quantum physics, 3rd ed. (John Wiley & Sons, 2003). 

[58] R.W. Boyd, Nonlinear Optics, 3rd ed. (Academic Press, 2008). 

[59] J.D. Cox and F.J. Garcia de Abajo, Nat. Commun. 5, 5725 (2014). 


